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Abstract 



In this paper we continue and extend a systematic study of plateaux in magne- 
tization curves of antiferromagnetic Heisenberg spin-1/2 ladders. We first review 
a bosonic field-theoretical formulation of a single XXZ-chain in the presence of a 
magnetic field, which is then used for an Abelian bosonization analysis of A'^ weakly 
coupled chains. Predictions for the universality classes of the phase transitions at 
the plateaux boundaries are obtained in addition to a quantization condition for 
the value of the magnetization on a plateau. These results are complemented by 
and checked against strong-coupling expansions. Finally, we analyze the strong- 
coupling effective Hamiltonian for an odd number of cylindrically coupled chains 
numerically. For = 3 we explicitly observe a spin-gap with a massive spinon-type 
fundamental excitation and obtain indications that this gap probably survives the 
limit N ^ oo. 
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I. INTRODUCTION 



So-called "spin ladders" have recently attracted a considerable amount of attention (for 
reviews see e.g. [1]). They consist of coupled one-dimensional chains and may be regarded 
as interpolating truly one- and two-dimensional systems. Such an intermediate situation 
may be useful (among others) for the understanding of high-Tc superconductors. In fact, 
modifications of the high- Tc materials (see e.g. [2]) give rise to experimental realizations of 
spin ladders. However, the field was motivated by an observation that mainly concerns the 
magnetic spin degrees of freedom, namely the appearance of a spin gap in N — 2 coupled 
gapless chains (see e.g. [3]). 

One-dimensional quantum magnets have been studied in great detail over the past 
decades. One remarkable observation in this area is the so-called "Haldane conjecture" 
[4] which states that isotropic half-integer spin Heisenberg chains are gapless while those 
with integer spin are gapped. Although this statement has not been proven rigorously yet, 
a wealth of evidence supporting this conjecture has accumulated in the meantime [5] (see 
also [6] for a recent field theoretical treatment). 

Spin ladders are more general quasi one-dimensional quantum magnets. Again, one of 
the attractions is a natural generalization of Haldane's conjecture [4] to such N coupled 
spin-S" chains: If SN is an integer, one expects a gap in zero field, otherwise not. This 
conjecture is suggested among others by the large- limit (see [7] for a recent review and 
references therein), strong-coupling considerations [3,8], numerical computations [9-11] and 
even experiments (see e.g. [2]) ^. If one includes a strong magnetic field, these Haldane 
gaps become just a special case of plateaux in magnetization curves. In the presence of a 
magnetic field, one of the central issues is the quantization condition on the magnetization 
(M) for the appearance of such plateaux, which for the purposes of the present paper will 
be some special form of 

ISN{1- (M)) eZ. (1.1) 

Here (M) is normalized to a saturation value (M) = ±1 and I is the number of lattice sites 
to which translational symmetry is either spontaneously or explicitly broken. Haldane's 
original conjecture [4] is related to I — N = 1, (M) = 0. More general cases ior N — 1 were 
treated in [12-14]. In an earlier paper [15], we have studied realizations of the condition 
(1.1) for N > 1 but with the specializations / = 1 and S" = 1/2. 

So far, spin ladders in strong magnetic fields have attracted surprisingly little attention: 
To our knowledge, only the case of an = 2 leg ladder had been investigated prior to [15]. 
The experimental measurement of the magnetization curve of the organic two-leg ladder ma- 
terial Cu2(C5Hi2N2)2Cl4 [16] gave rise to theoretical studies using numerical diagonahzation 
[17], scries expansions [18] and a bosonic field theory approach [19]. In this case, {i.e. N = 2 
and iS = 1/2) there is a spin gap which gives rise to an (M) = plateau in the magnetization 
curve. The transition between this zero magnetization plateaux and saturation is smooth 
and no non-trivial effects (in particular no symmetry breaking) were observed. 



This result was obtained for planar configurations and often even an equal strength of coupling 
constants is assumed. Such properties may be crucial as we will show below. 
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For N > 2 one can expect plateaux at non-trivial iV-dependent fractions of the saturation 
magnetization [15]. Though this was a new observation for spin ladders, the phenomenon 
itself is not completely new. For example, the possibility of magnetization plateaux for the 
case of single spin- 5" Heisenberg chains has been discussed systematically in [12] which also 
motivated some of our work. The attraction of this phenomenon in spin ladders is that they 
provide clear and natural realizations of such plateaux. For example a numerical analysis of 
the case = 3 explicitly exhibits a robust plateau with (M) = 1/3 [15] which should also 
be observable experimentally, e.g. in a suitable organic spin-ladder material. 

It is the purpose of the present paper to continue a systematic study of (1.1) for generic 
N > 1 hy using different complementary techniques, such as Abelian bosonization, strong 
coupling expansions, and numerical computations (the reader may find e.g. [15] helpful in 
the understanding of the present work). Here, we will among others provide evidence that 
for = 1/2 spontaneous breaking of the translational symmetry to / = 2 can be induced 
by strong frustration or an Ising-like anisotropy, while I > 3 presumably needs explicit 
symmetry breaking. (In [14] a slightly different example of spontaneous symmetry breaking 
with I = 2 was studied). 

The prefactor ISN in (1.1) may seem quite cumbersome, but it just counts the possible 
5'^-values in a unit cell of a one-dimensional translationally invariant groundstate. It should 
also be noted that (1.1) is just a necessary condition; whether a plateau actually appears 
or not depends on the parameters and the details of the model under consideration. For 
example, plateaux with non-zero (M) ^ have not been observed in the SU{2) symmetric 
higher spin-S* Heisenberg chains (see e.g. [20]), unless translational invariance is explicitly 
broken (c.f. [13] for S = 1). 

Conditions of the type (1.1) occur also in generalizations of the Lieb-Schultz-Mattis 
theorem [21-23,12]. This theorem constructs a non-magnetic excitation which in the ther- 
modynamic limit is degenerate with the groundstate for a given magnetization (M) and 
orthogonal to it unless (M) satisfies (1.1). In this manner one proves the existence of either 
gapless excitations or spontaneous breaking of translational symmetry. Unfortunately it is 
at present not clear that this theorem applies to plateaux in magnetization curves since they 
require a gap to magnetic excitations. 

Here we concentrate mainly on the case S — 1/2 and all couplings in the antiferromag- 
netic regimes, but try to keep as general as possible. Other situations can be analyzed 
as well, but may lead to somewhat different physics (compare e.g. [24] for the example of 

= 3 antiferromagnetically coupled ferromagnetic chains). 

This paper is organized as follows: In section II we first review some aspects of the 
formulation of a single XXZ-chain as a bosonic c — 1 conformal field theory. This serves 
as a basis of later investigations and illustrates some generic features also present in A^-leg 
spin ladders. In section 111 we first introduce the precise lattice model and its field-theoretic 
counterpart which we then analyze in the weak-coupling regime. Section IV starts from 
the other extreme -the strong-coupling limit- and proceeds with series expansions around 
this limit. In section V we discuss an effective Hamiltonian for the strong-coupling limit 
of an odd number N of cylindrically coupled chains which we then analyze numerically in 
section VI. We summarize our results by presenting "magnetic phase diagrams" in section 
VII before we conclude with some comments and open problems (section VIII). 
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II. A SINGLE XXZ-CHAIN 



First, we recall some results for the XXZ-chain on a ring of L sites in the presence of a 
magnetic field h applied along the z-axis: 

Hxxz ^ ^S^S^^i + 9 {S^Sx+i + 'S'^'S'j+i) \ -h^S^. (2.1) 

x=l ^ ^ x=l 

Apart from being the basis for the investigation in later sections, this also serves as an 
illustration of some general features. It should be noted that in (2.1) the magnetic field 
is coupled to a conserved quantity which is related to the magnetization (M) via (M) = 
J J2x=i ^xY I*'^^ t^is reason, properties of (2.1) in the presence of a magnetic field h 
can be related to those a,t h — and the magnetic field h can be considered as a chemical 
potential. 

The Hamiltonian (2.1) is exactly solvable by Bethe ansatz also for h ^ 0. In this way it 
can be rigorously shown that its low-energy properties are described by a c = 1 conformal 
field theory of a free bosonic field compactified at radius R in the thermodynamic limit 
for A > —1 and any given magnetization (M) (see e.g. [25] and compare also [26] for a 
detailed discussion of the case |A| < 1). More precisely, upon insertion of the bosonized 
representation of the spin operators into the Hamiltonian (2.1) (see e.g. [27]) one obtains 
the following low-energy effective Hamiltonian for the XXZ-chain: 

Hxxz = / dx^ {U\x) + R\{M),A) (9.0(x))^} (2.2) 

with H = ^dx4>, and (p = (pL + (pR, (p = (pL — <Pr- In (2.2) we have suppressed a for our 
purposes irrelevant proportionality constant that includes the velocity of sound. In this 
formulation, the effect of both the magnetic field h and XXZ-anisotropy A turns up only via 
the radius of compactification R{{M),A). This radius governs the conformal dimensions, 

in particular the conformal dimension of a vertex operator e'^*^ is given by (4:^) ■ We now 
describe how R can be computed. 

We parametrize the XXZ-anisotropy by A = cos^ with < ^ < tt for — 1 < A < 1 
and by A = cosh 7 with 7 > for A > 1. Now for given magnetization (M) > and 
XXZ-anisotropy A, the associated radius of compactification R and magnetic field h/ J can 
be obtained by solving integral equations (see e.g. [26,28-30]) in the following way: Firstly, 
introduce a function a{r]) for the density of particles satisfying the integral equation 

(t{v) = ^ l^giv) - J^^ K{7j - v'Hv')dv''^ (2.3) 

where the kernel K(rj) and the righthand side g{rj) are presented in Table I. 

The real parameter A > in (2.3) describes the spectral parameter value at the Fermi 
surface and is determined by the magnetization (M) via the filling condition 

|\(r^)dr^ = i(l-(M)) . (2.4) 

In general, one has to adjust A iteratively by first numerically solving (2.3) and then checking 
for (2.4). Only some special cases can be solved exphcitly. This includes the case (M) — 
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A 


K{v) 






eo(^) 


cos (9 = A < 1 

A = 1 
cosh 7 = A > 1 


tan 6 


cot 1 


ft 


sin^ 


tan2 6» cosh^ a+sinh^ 1 
4 

??^+4 
tanh7 


cosh2 H+cot2 1 sinh2 2 
2 

coth ^ 


J 
ft 


cosh r;— cos 9 
ft 2 

sinh^ 7 


tanh"^ 7 cos^ ^+s'm'^ 2 


cos2 l+coth^ ^ sin^ 2 


J 


cos rj— cosh 7 



TABLE I. Functions appearing in the integral equations. 



and A < 1 where A = cxo is the correct choice. Once the desired value of A is determined, 
one introduces a dressed charge function ^(77) (see e.g. [28,29]) as a solution of the integral 
equation 

giving directly rise to the radius of compactification 

RUM), A) = 

If one further wants to determine the associated magnetic field h, one has to introduce 
another function ed{r]) , the dressed energy, satisfying the integral equation 

ed{r]) = eoiv) - ^ J^^ K{rj - v')ed{v')dv' (2.7) 

with the bare energy eo(r^) hsted in Table I. Then the magnetic field h/J is determined by 
the condition that the energy of the dressed excitations vanishes at the Fermi surface 

6d(A) = 0. (2.8) 

Using (2.5) one can sec that e{r]) = e{r])\h=o + jC{v) solves (2.7) if €{r])\h=o solves (2.7) with 
formally h = (but for the given (M)). From this and (2.8) one can easily obtain h once A 
is known via h/ J = — e(A)|/J=o/C(■A)• 
In general, these integral equations have to be solved numerically and A has to be de- 
termined by some iterative method. Although this is readily done by standard methods, a 
generally accessible implementation seems to be still unavailable. We have therefore decided 
to tentatively provide access to such solutions on the WWW [31]. This implementation 
works in the way described above. Typically, it gives results with an absolute accuracy of 
10~^ or better. Of course, one can change the order of the procedure: For example one could 
also prescribe h/J, then determine A from (2.7) and (2.8), next the radius of compactifica- 
tion R from (2.5) and (2.6) and finally (optionally) the magnetization (M) from (2.3) and 
(2.4). 

Some remarks are in order concerning the case A > 1: In this region the identification 
with a c = 1 conformal field theory has been established directly only by a numerical 
analysis of the Bethe-ansatz equations (for a summary see e.g. section 2 of [32]). However, 
the six-vertex model in external fields has been shown [33] to yield a c = 1 conformal field 
theory, and since the XXZ-chain arises as the Hamiltonian limit of the six-vertex model, this 
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(2.5) 



(2.6) 



indirectly establishes the identification. Still, explicit formulae e.g. for R are not available in 
the literature and therefore we have obtained the integral equations presented for the case 
A > 1 above by an analytical continuation of those for A < 1, as is suggested e.g. by [34] 
(see also [35] - apparently such a continuation was also used in the recent work [13]). We 
have performed some checks that this yields indeed correct results. For example, some radii 
R associated to certain magnetic fields h and anisotropics A obtained numerically for chains 
of length up to L = 234 [36] are reproduced in this way. 

As a further check, one can compare the critical magnetic field for the boundary of the 
(M) = plateau at A > 1 obtained by numerical solution of the above integral equations 
with the exact solution of the Bethe-ansatz equations [37] 

he _ 27rsinh7 ^ 1 

T ~ ^ 2^ , (2n+l)7r2 ' V^-^^ 

•J 7 n=o cosh 

where as before A = cosh 7. In this case (i.e. for (M) = 0) one has A — n and the above 
integral equations can be solved using Fourier series [34,35,37]. If one solves the above 
integral equations numerically, the deviation from the exact result (2.9) is of the order of 
the numerical accuracy (in our implementation [31] always less than 10~^). 

In passing we make a comment which will turn out to be useful later. The result (2.9) 
is the gap to 5"^ = 1 excitations. However, the fundamental excitation of the XXZ-chain 
is known to be a so-called "spinon" which carries — 1/2 [38-40]. This spinon can be 
regarded as a domain- wall between the two antiferromagnetic groundstates for A > 1. Since 
a single spin-fiip creates two domain-walls, the lowest = 1 excitation is a scattering state 
of two spinons. This picture can be useful e.g. in numerical computations. For example, a 
single spinon can be observed for odd L with periodic boundary conditions. 

After this digression let us now return to the above integral equations. The results 
obtained from them are summarized in the magnetic phase diagram for the XXZ-chain Fig. 
1 (see also [41,32] for similar pictures). There are two gapped phases: A ferromagnetic one at 
sufficiently strong fields (which is actually the only one for A < —1) and an antiferromagnetic 
phase for A > 1 at small fields. In between is the massless phase where the bosonized form 
(2.2) is valid. An elementary computation of the spinwave dispersion above the ferromagnetic 
groundstate shows that the transition between the ferromagnetic phase and the massless 
phase is located a.t h/ J = 1 + A. This transition is a very clear example of the DN-PT 
universality class [42,43] ^, i.e. for (M) — 1 the magnetization behaves as (compare also 
[44]) 

{{M) - Mef ^ - hi (2.10) 
with here Mf. — 1 and hj J = 1 -|- A. At this transition the radius takes the universal value 



^This terminology is motivated as follows: This type of transition is usually known under the name 
"Pokrovsky-Talapov" [43], but in the context of magnetization processes was actually discovered 
first by Dzhaparidze and Nersesyan [42]. 
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FIG. 1. Magnetic phase diagram of the XXZ-chain (2.1). For explanations compare the text. 

The other transition hne starts at the SU{2) symmetric point (M) = 0, A = 1 with a 
radius i?(0, 1) = Actually, at (M) = the additional operator 

cos (4V7r0) (2.11) 

appears in the continuum limit, which we have suppressed in (2.2) since it is irrelevant 
inside the massless phase. At A = 1 it is marginal and becomes relevant for A > 1, opening 
the gap that gives the boundary of the antiferromagnetic phase in Fig. 1. The associated 
phase transition is a Kosterlitz-Thouless (K-T) transition [45] (see e.g. [46-48]). The almost 
marginal nature of the operator responsible for the gap leads to a stretched exponential 
decay for A shghtly bigger than one which is characteristic for a K-T transition [49]. The 
exact asymptotic form for the gap (or critical magnetic field) of the XXZ-chain is easily 
obtained from the Bethe-ansatz solution (2.9) upon noting that 7 ~ y2(A — 1) and that in 
the limit A ^ 1 only the term n — contributes to the sum. One then finds [37]: 

h .-^ 

_^ ^ 47re 2V^(A^ (for A slightly bigger than 1). (2.12) 

For this reason the phase boundary is indistinguishable from the h — Q line for XXZ- 
anisotropies up to A ^ 1.2 on the scale of Fig. 1. In this region, the numerical determination 
of the radius R is difficult for (M) — > 0. Nevertheless, using that A = tt for (M) = and 
A > 1, one can readily check that the constant function ^[r]) = \ solves the integral equation 
(2.5). Then one obtains from (2.6) that i?(0, A > 1) = ^. 

For (M) = and |A| < 1 one has A = cxo (see above) and one can use the Wiener- 
Hopf method to solve the above integral equations in closed form. This yields in particular 
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R{0, A) = (l - ^ cos-i Aj [27]. In general, the radius R increases with increasing A. 

For A > 0, it decreases with increasing magnetization (M), while for A < this is reversed 
to an increase in R with increasing (M) . Clearly, the radius must be constant on the XY-hne 
which separates these two regions: R{{M),0) = (though the magnetic field associated 
to a given (M) is still a non-trivial function which should be computed from the above 
integral equations). 



<M> 




FIG. 2. Magnetization curve of the XXZ-chain at A = 2 obtained from the integral equations 

(2.3), (2.4), (2.7) and (2.8). The inset shows the region of small magnetization and illustrates that 
also the transition (M) ^ is compatible with the DN-PT universality class for A > 1 (the dashed 
line is a fit to the universal form (2.10)). 

Including the operator (2.11) in the bosonized language (2.2), one recovers a Hamiltonian 
treated in [50] as a model for commensurate-incommensurate transitions. This means that 
the transition (M) — for A > 1 is predicted to be in the DN-PT universality class 
[42,43], too. The same bosonization argument also leads to the already mentioned result 
i?(0,A > 1) = (see also [25]). The inset in Fig. 2 illustrates for A = 2 that for 
sufficiently small magnetizations one can indeed observe a behaviour that is compatible 
with the universal square- root (2.10). It should be noted though that the window for the 
universal DN-PT behaviour is too small to permit verification within our numerical accuracy 
for smaller values of A (e.g. A < 1.2) where neither a reliable numerical check of the result 
R — ^ is possible. An analytic check of the asymptotic behaviour of the magnetization 



8 



from the Bethe-ansatz solution would be interesting, but is beyond the scope of the present 
paper. 

It should be noted that the height of the entire inset in Fig. 2 corresponds to the first 
step in the magnetization curve of a chain of the finite size L — 112. Therefore, the exact 
solution is crucial in verifying the exponent (2.10) - a numerical or experimental verification 
of this behaviour restricted to such a small region would be extremely difficult. Similar 
difficulties will be faced in an experimental or numerical verification of -R(0, A) = or the 
equivalent statement for the correlation function exponents in the region of A slightly larger 
than 1. 

III. LADDERS: WEAK COUPLING AND ABELIAN BOSONIZATION 

In this section we apply Abelian bosonization to the weak-coupling region J' <^ J of 
iV-leg spin ladders. In particular, we will show how the necessary condition (1.1) arises 
in this formulation and discTiss Tinder which circumstances an allowed plateau does indeed 
open as a function of the parameters J' and A. The lattice Hamiltonian for this system is 
given by 

i,j x=l 1=1 x=l i,x 

(3.1) 

where J' and J are respectively the interchain and intrachain couplings, h is the external 
magnetic field and the indices i and j label the different chains (legs) in the ladder. The 
sum in the first term is over all possible couplings between chains. The case of periodic 
boundary conditions (PBC) and open boundary conditions (OBC) will be discussed later. 
Here we have explicitly included an XXZ-anisotropy A in the intrachain coupling. We have 
kept the interchain coupling J' SU{2) symmetric for simplicity in later sections although 
this is not substantial in the weak-coupling region which we will discuss in the remainder of 
this section. 

The corresponding effective field-theoretic Hamiltonian is obtained using standard meth- 
ods [51,27] (see also [12-14] for the case of non-zero magnetization). One essentially uses 
(2.2) as the effective Hamiltonian for each chain and the bosonized expressions for the spin 
operators which read: 

^ + canst. : cos(2A;^x + ^4^0^) : , (3.2) 

and 

St,x ~ : e^^^^*(l + const. cos{21^pX + ViTr^O) : ■ (3.3) 

Here we have set a lattice constant to unity which appears in passing to the continuum 
limit. The colons denote normal ordering which we take with respect to the groundstate 
of a given mean magnetization (Mj) in the ith chain which is a natural choice. This leads 
to the constant term in (3.2) which will play an important role in the discussion of the 
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= / da; 



terms that can be generated radiatively. The prefactor 1/2 arises from our normahzation 
of the magnetization to saturation values (M) = ±1. The Fermi momenta kp are given by 
ki, = 7r(l - {Mi))/2. 

In the weak-coupling limit along the rungs, J' <C J, we obtain the following bosonized 
low-energy effective Hamiltonian for the A^-leg ladder keeping only the most relevant per- 
turbation terms: 

- N 

+ ^|a2 : cos(2x(A;i. + 4) + ^4^(0^ + 0,)) : (3.4) 

-FAs : cos(2a;(/c^ - k{r) -\/47r(0i - 0^)) : : cos (\/n{4>i - ^ 

The four couphng constants Aj essentially correspond to the coupling J' between the chains: 
Ai ~ J' /J. In arriving to the Hamiltonian (3.4) we have discarded a constant term and 
absorbed a term linear in the derivatives of the free bosons into a redefinition of the applied 
magnetic field. 

The Hamiltonian (3.4) has been also used to represent spin-A^/2 chains (see e.g. [51,52]), 
since they can be obtained in the limit of N strongly ferromagnetically coupled chains 
(J' — s> — oo). However, here we will analyze (3.4) mainly in the case of small antiferromagntic 
J' and discuss various boundary conditions. 

Note that the A2 and A3 perturbation terms contain an explicit dependence on the posi- 
tion (in the latter case this x-dependence disappears for symmetric configurations with equal 
kp). Such operators survive in passing from the lattice to the continuum model, assuming 
that the fields vary slowly, only when they are commensurate. In particular, the A2 term 
appears in the continuum limit only if the oscillating factor exp{i2x{kp + kp)) equals unity. 
If the configuration is symmetric, this in turn happens only for zero magnetization (apart 
from the trivial case of saturation) . 

For simplicity let us first analyze the case with N = 3 and PBC We first have to 
diagonalize the Gaussian (derivative) part of the Hamiltonian. This is achieved by the 
following change of variables in the fields: 

11 1 

V'l = ^ (01 - 03) , = ^ (01 + 03 - 202) , -00 = (01 + 02 + 03) ■ (3.5) 

In terms of these fields the derivative part of the Hamiltonian can be written as: 

^der. = / dx| {r\{M),A) [(1 + a) {d^iPDix))' + {i-b) {{d,iPi{x))^ + {d.Mx))')] } 

(3.6) 

where a = 2J' / {Jrc^R"^) = 2b. We can now study the large-scale behaviour of the effective 
Hamiltonian (3.4) where we assume all kp equal due to the symmetry of the chosen config- 
uration of couplings. Let us first consider the case when the magnetization (M) is non-zero. 
In this case only the A3 and A4 terms are present. The one-loop R.G. equations are: 
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dlnL \ 2B? 

^ = (2-2rf(l-6))A4-7^A^ (3.7) 

It is important to notice that only the fields V^i and %Ij2 enter in these R.G. equations, since 
the perturbing operators do not contain the field The behaviour of these R.G. equations 
depends on the value of R. The main point is that always one of the two A perturbation 
terms will dominate and the corresponding cosine operator tends to order the associated 
fields. This gives a finite correlation length in correlation functions containing the fields 
ipi and %p2 (or their duals). For example, for A < 1 we have that < (27r)~^ since 
(M) 7^ 0. Then, from (3.7) one can easily see that the dominant term will be the A4 one. 
This term orders the dual fields associated with V'l and ■02- Then, the correlation functions 
involving these last fields decay exponentially to zero. In either case, the field ipD remains 
massless. A more careful analysis of the original Hamiltonian shows that this diagonal field 
will be coupled to the massive ones only through very irrelevant operators giving rise to 
a renormalization of its compactification radius. However, due to the strong irrelevance of 
such coupling terms these corrections to the radius are expected to be small, implying that 
the value of the large-scale effective radius keeps close to the zero-loop result R^/l — a. It 
is straightforward to generalize this to chains when all possible coupling are present and 
have the same value J'. One can find a change of variables on the fields to: 

ipD ; ipi i^l,...,N-l 

where ^^d — l/V^Z^ili 0i- Again, for non-zero magnetization, all but the diagonal field 
■0D will be present in the perturbing terms A3 and A4. The R.G. equations are essentially 
the same as (3.7) and the result is that only the field ipD will be massless. 

Wc arc then left in principle with a free Gaussian action for the diagonal field. However 
some operators can be radiatively generated. We see from eqs. (3.2,3.3) that when we turn 
on the interchain coupling, the "N-Umklapp" term 

2x^4 + v^^(/)J = J'^ coshx^^kp + y^^M (3.8) 

i=l i=l J V 1=1 / 

appears in the operator product expansion (OPE). 

Again, this operator survives in passing from the lattice to the continuum model, assum- 
ing that the fields vary slowly, only when the oscillating factor exTp{i2xJ2j=i kp) equals one. 
This in turn will happen when the following specialized version of the condition (1.1) 

^(l-(M))eZ (3.9) 

is satisfied. At such values of the magnetization, the field ipD can then undergo a K-T 
transition to a massive phase, indicating the presence of a plateau in the magnetization 
curve. An estimate of the value of J' at which this operator becomes relevant can be 
obtained from its scaling dimension which in zero-loop approximation is given by 
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dim (cos (74^^.)) = ^^^^^^^ . (3.10) 

At A = 1 one then obtains ^ 0.09J for the (M) = 1/3 plateau at iV = 3 and ^ 0.7J 
for (M) = l/2atA'" = 4 and also for (M) = 1/5 at = 5. At the opening of such plateaux, 
the effective radius of compactification is fixed to be: 

^eff. = £ (3.11) 

and the large-scale effective spin operators are (c.f. [51]): 

'S'eff {x) fti J + const. : cos(2A;j.a; + \fA^^PD) ■ +^ , (3.12) 
V 27r ox 2 

and 

'S'eff.(^) ~ : e^iV^^o^i + coj^gt cos{2kFX + Vi^ipo)) ■■ ■ (3.13) 
Then, eq. (3.11) fixes the values of the correlation exponents at this point to be 

?7z = 4; Vxy^^- (3.14) 

On the other hand, commensurate-incommensurate transition results [50,19,13] imply that 
the values of these exponents should be 

^z = 2; Vxy^^ (3.15) 

along the upper and lower boundary of a plateau. This situation is similar to the XXZ-chain 
at A = 1 and A > 1 for the boundary of the (M) = plateau. However, the values of the 
exponents are different since the perturbing operators are different. 

Note that the "N-Umklapp" process which allows the appearance of (3.8) produces a 
complete family of operators given by: 

/ iV N \ / N \ 

COS 2x1 J2ki, + = cos 2x1 Y.^f + IV^ttNt/jd (3.16) 

\ 1=1 i=l J \ i=l / 

with / an arbitrary integer. The values of the magnetization for which one of these operators 
is allowed are subject to a generalization of (3.9), namely (1.1) in the Introduction (with 
S — 1/2). However, the dimensions of these operators increase with P. So, these operators 
cannot be relevant unless we consider regimes with an anisotropy parameter A bigger than 
one or very big values of the interchain coupling J' far from the perturbative regime of the 
present analysis. Therefore, higher values of / are realized only under special conditions. 
While / = 2 can be obtained by either strong Ising-like anisotropy A or frustration at strong 
coupling (see section VI below), it is possible that / > 3 can be realized only if suitable 
symmetry breaking terms are explicitly introduced into the Hamiltonian (3.1). 
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Note that formally, the preceding analysis can also be carried out using the fermionic 
Jordan- Wigner formulation. For example, in this formulation the "N-Umklapp" operator 
(3.8) is given by 

where Ra and La are the right- and left-moving components of the fermions. We have 
chosen to use the bosonized language because it is more appropriate for general values of 
the anisotropy A. 

The analysis above was for the case where all the chains were coupled together with the 
same coupling value. More precisely, the estimates for the appearance of plateaux were for 
positive (frustrating) interchain couphng. To generalize this to PBC (which is different from 
the preceding case for N > 4). wc first notice, using the bosonized expression of the effective 
Hamiltonian, that this configuration of couplings is not stable under R.G. transformation. 
E.g. the OPE between terms like cos(0i — ^2) and cos(02 — ^3) generates an effective coupling 
between the fields (pi and 03, etc. The underlying intuitive picture is that antiferromagnetic 
couplings between the chain 2 with the chains 1 and 3 generates an effective ferromagnetic 
coupling between the chains 1 and 3. For example, for = 4 and PBC, ferromagnetic 
couplings are generated along the diagonals between originally uncoupled chains. This 
case is part of the family of configurations with antiferromagnetic nearest neighbour and 
ferromagnetic next-nearest neighbour couplings. For this general situation at A'^ = 4, the 
coupling matrix in the derivative part is given by: 



(3.17) 



where a and b are positive. As in the preceding analysis, one can change variables to 

'ipD = + 02 + 03 + 04) ; A = ;^(0i - 02 + 03 - 04) ; 

V'2 = ^(01 - 03) ; ^3 = ^(02 - 04). (3.18) 

For generic values of the magnetization, it is easy to see that the diagonal field is again 
the only field that does not acquire a mass under the perturbation. Then, the analysis of the 
appearance of the N-Umklapp term for particular values of the magnetization is identical 
to the one performed before. The generalization to generic N with PBC is straightforward, 
one first builds the radiatively generated couplings by keeping only the lowest order in J'. 
Once this step is performed, the only difference with respect to the case of equal interchain 
couplings is the zero-loop value of the dimension of the N-Umklapp operator (which enters 
via the initial conditions for the R.G. fiow). This has the effect of changing the value of 
the coupling J' at which a plateau opens with a given value of the magnetization, but the 
qualitative behaviour of the system is similar. This conclusion is not so straightforward 
for (M) = 0, where as we will see, the difference between frustrating and non-frustrating 
configurations can become crucial. 



/ 


1 


a 


-b 


a \ 




a 


1 


a 


-b 




-b 


a 


1 


a 


V 


a 


-b 


a 


1 / 
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Concerning finally the case of OBC, let us first consider again the case = 3 with 
antiferromagnetic coupling between the first and second chain and the second and the third 
chain. Again, this coupling is not stable under R.G. transformation. Under R.G. trans- 
formations the OBC configuration will flow towards a non-frustrating cyclically coupled 
configuration. The main point is that for weak coupling and non-zero magnetization, the 
most relevant perturbing term will be again the one containing differences of fields or their 
duals. Then they will produce a mass gap for all the relative degrees of freedom and one 
recovers a scenario similar to the symmetric case, where only one massless field is left. On 
the other hand, the appearance of "N-Umklapp" operators and their commensurability is 
unchanged, since these criteria depend on the value of the magnetization and not on the 
particular couplings between the chains. 

Let us study now the more complicated case of zero magnetization. For (M) = the 
A2 term in (3.4) is commensurate and must be included in the perturbation terms. The 
situation is now much more complicated because this relevant operator couples the diagonal 
field ■0D with the massive ones. For equal coupling between N — 3 chains, the R.G. equations 
are now: 



da ISttXI 



dlnL 

db . ( \l 3A? 



d\nL 



dXo /„ 2 1 / 2 1 



dlnL 
dXs 



dlnL 

dX 



dlnL 

with the R.G. initial conditions 




+ Un^R^X, 



4 



Ao — TtAoA; 



2^3 



A3 — 7rA3 — 7rA2 



TrXi (3.19) 



2 J' 

a(0) = 26(0) = (3.20) 

and 

A2(0) = A3(0) = l/2{constyj'/J ; A4(0) = J'/ J (3.21) 

where we kept the notation of (3.4,3.6). We see that the radius of compactification of the 
diagonal field is now strongly affected by the presence of the A2 term. Note also that the 
"N-Umklapp" process generates the operator 

J'^ cos (yA^i/jo) (3.22) 

for N even, and 

J'^ cos (2\/4^^ip) (3.23) 

for N odd. For non-frustrating interchain coupling (a negative J' coupling between all the 
chains for example), all relative fields are massive according to [52]. We can then integrate 
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out these massive degrees of freedom. The crucial point is that now the radius of the diagonal 
field gets a non-trivial correction due to the strong interaction with the massive fields. Since 
this field is the only one expected to describe the large-scale behaviour of the system, for 
A = 1 and (M) = 0, the SU (2) symmetry of the model would fix the radius of this field to 



For such a value of the renormalized radius, the "N-Umklapp" term becomes strongly rel- 
evant for even, and marginally irrelevant for odd. These arguments are based on the 
assumption that the (uncontrolled) R.G. flow will drive our system to the {SU (2) symmetric) 
strong coupling regime. The situation is even more subtle for positive J' (or Aj), because in 
this case, from (3.19) one sees that the quadratic terms could now prevent the R.G. flow to 
reach the same strong coupling regime as for J' < 0. A numerical analysis of the R.G. flow 
for a frustrated three-leg Hubbard ladder at half filling provides evidence for the opening of a 
gap [53]. On the other hand, a non-abelian bosonization analysis [54] lead to the conclusion 
that the weak-coupling region is gapless. This case deserves further investigation and series 
expansions are one way to approach this issue. 

IV. STRONG COUPLING EXPANSIONS FOR AT-LEG LADDERS 

In this section we diagonalize the interaction along the rungs exactly for J = and then 
expand quantities of interest in powers of J/ J' around this limit. In order to be able to cover 
a variety of cases, we used a quite general method to perform the series expansions which 
is summarized e.g. in section 3 of [55] (actually, the program used in the present paper is a 
modifled version of the one used loc.cit.). 

As was already pointed out in [15], one can simply count the number of chains N in the 
limit J/J'-^O in order to determine the allowed values of the magnetization (M). This is 
presumably the simplest way to obtain the quantization condition (3.9). A less trivial fact is 
that all these values of the magnetization are in fact realized. For example, for ferromagnetic 
coupling J' < 0, the magnetization jumps immediately from one saturated value (M) = — 1 
to the other one ((M) = +1) as the magnetic field h passes through zero. Nevertheless, 
for not too large N one can readily compute the magnetization curve of (2.1) and check for 
antiferromagnetic coupling J' > that all possible values of the magnetization are indeed 
successively reahzed as the field is increased. The critical magnetic fields h at which one 
value of the magnetization jumps to the next largest one arc given in Table II. 

As a next step, one can take the intrachain coupling J perturbatively into account. 
First, we look at a two-leg ladder {N = 2). The rung Hamiltonian Hr = J' 5*1 5*2 has two 
eigenvalues whose difference corresponds to the critical fields presented in Table II. The lower 
eigenvalue equals — 3J'/4 and belongs to the spin S — eigenstate while the other threefold 
degenerate one equals J'/4 and corresponds to the spin triplet {S = 1). For convenience, we 
concentrate on an isotropic interaction for the rungs, but it is straightforward to include an 
XXZ-anisotropy A in the interaction along the chains. One motivation for doing so is that 
this permits further comparison with the weak-coupling analysis (J' <^ J) of the previous 



be [52]: 




(3.24) 
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2 
3 
4 
5 
6 


±-, 

= ±1.80902, ±1.11887, 
±1.86603, ±1.38597, ±0.49158 


±- 
±2, ±1 

_l_5+\/5 _i_3+V5 n 
-■- 4 ' -■- 4 ^ ^ 

_i_o _l1+\/5 i\/T3-\/5 
=nz, in 2 ) =■= 2 




TABLE II. Values of h at which the magnetization jumps for (2.1) 


with coupling constant J', 



A = 1, iV sites and different boundary conditions. 



section. At J = the groundstate is obtained by putting singlets on each rung. A basic 
excitation at J = is given by one triplet in a sea of singlets. Since the SU{2) symmetry is 
broken down to U{1) by the perturbation, different series are obtained for the 5"^ = ±1 and 
— components of the triplet. 

Here we concentrate just on the series for the gap, but also previous results for the 
groundstate energy and the dispersion relations are readily extended to higher orders [8] , or 
to analytical expressions in A for longer series [56] at A = 1 with numerical coefficients. 

The gap is obtained by the value of the excitation energy of a single flipped spin at 
momentum k = n with = ±1. We find 

E2 , fJ\ 1 + AV^\^ (l + AffJ\^ -2 + 6A-9A2 + AV^ 

21 - 84A + 39A2 - 48A3 + 2A^ f 
H 



256 VJ'. 
82 - 98A + 155A2 - 50A=^ + 80A^ - 12A^ / J ^ ^ 




+ . (4.1) 



1024 \J' 

At the isotropic point A = 1 we recover well-known results: For this special case, the first 
three orders can be found in [8], a fourth order was given in [15] and numerical values of the 
coefficients until 13th order are contained in [56]. 

The series (4.1) contains a singularity at J' = which has no physical meaning but is 
simply due to the choice of expansion parameter. We therefore analyze it by removing this 
singularity via the substitutions 

-=7TJ • * = '-'(j7)- (^-2) 

Prom the raw transformed series one can then find some indication of an extended massless 
phase at small J' if A < Ac with Ac ~ 0.25 . . .0.5. The opening of this massless phase is 
predicted by the zero-loop analysis of the previous section to take place at Ac = 0. Since 
the information obtained in the weak-coupling regime from a strong-coupling series is not 
extremely accurate, this agreement can be considered reasonable. 

Now we turn to = 3 and OBC. In a way similar to the previously discussed series one 
finds the following fourth order series for the lower and upper boundary of the (M) = 1/3 
plateau: 
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^ ^ fA . n Z _ (A + l) (8A-5) fJ_\^ (A + l)(142A^- 307A-23) / J_\^ 
J' ^ ' J' 27 VJV 972 [j'J 

(A+l)(40572A3- 83025A2 + 7696lA- 73295) / .^(fJ\^\ 
^ 367416 [jj ^ J ' ^ ' 

/ic2 _ 3 J 10 + 17A2/J\2 2196A + 252 - 554A3 + 171A2 / J\3 

7^ " 2 ~ 7 ^ 36 [l') ^ 3888 V 7 j 

30172 + 38988A - 28387 A^ + 7028A'^ - 8886A^ / J \ ^ o ( f 
^ 326592 [j') ^ [[j) ) ' ^ ■ 

A third-order version of these series was aheady presented in [15] for the special case A = 1. 
We employ again the transformations (4.2) to analyze these series. The raw transformed 
series indicate for A = 1 that the (M) = 1/3 plateau does not extend down until J' — but 
ends at a critical value J^. The numerical value is found to be ~ 1.0 . . . 1.4J at A = 1. 
This number should however not be taken too seriously as is also indicated by the large 
uncertainty of the critical anisotropy A^, above which this plateau extends over all non-zero 
J': Ac ~ 1.0 .. . 1.6. At least, this rough estimate for Ac is compatible with Ac 1.19 as 
obtained from the zero-loop weak-coupling analysis. 

The next case we consider is = 4 and PBC. In the strong-coupling limit we find 
plateaux at (M) = and at (M) = 1/2. Series can be computed readily for the gap (which 
determines the boundary of the (M) = plateau) and the lower and upper boundary of the 
(M) = 1/2 plateau. In this order, they read 

E'f^ ^ 4 /J\ 33A2-12A + 20 / J\^^ 24 + 194A2 + 131A /J\2 



J' 3 \J'J 108 VJV 1296 \J' 

^3524213 - 17599776A2 + 9014208 A + 1923768A3 + 7733988A4 / J ^ ^ 



39191040 \J' 




(4.5) 



/i^) , , 3A + 8 / J \ , 9 A^ + 96A - 308 / J 



J' 6 VJV 864 VJ' 

369A + 972A3 - 1314A2 - 9464 / J \ ^ 



2 



^ 31104 VJ'. 

885195A^ - 69076728A2 - 61318885 - 545832A3 + 117897360A / J ^ ^ 



156764160 VJ 




(4.6) 



^ , A-2 / J\ , 5A2 + 22 / J\^ 8A=^ - 42A - 9A2 - 51 / J 



2 + — ^K: + 



3 



J' 2 VJV 32 VJV 256 VJ'. 

^ 38A^ - 1981 A^ - 56A3 + 1634A - 403 / J y ^ c(('^)^^ (4 7) 



4096 \J'J \\J' 

The superscript '(p)' means that these series are for PBC. 

Again, we analyze these series using the transformations (4.2). We apply this first to 
the gap (4.5) and find that the gap closes for some J' > if A < Ac where the estimates 
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for the critical value span an interval A^. ~ 0.8 ... 1.2. This interval is centered around the 
value Ac = 1 predicted by power counting in the context of Abelian bosonization. 

Concerning the opening of the (M) = 1/2 plateau, we can first locate its ending-point 
in the same way as before as 0.8 . . . 1.6J at A = 1. What is more interesting is the 
conclusion that this ending-point cannot be pushed down to = by increasing A. This 
is in agreement with the zero-loop weak-coupling analysis which implies that an (M) — 1/2 
plateau does not exist for J' ^ J and = 4 regardless of the choice of A. 

Finally we present second order versions of analogous series for = 4 and OBC (denoted 
by a superscript 'o'): 



E 



J' 



^— f 1682A2 - 2760A + 1847) + ^— f 4176A2 - 7176A + 4063) 

1656 V / 3312 V ^ 



J' 



^ ^ - ^) - (^) - { li (---^ - - 155) 

f-1406A2 + 2691A + 637) + — — ( 27819A2 - 57408A - 97216) 
1656 V ) 26496 ^ ) 



+ 8k("^^^^^' + ^^^3A + 3920)|(j;)' + o(^(-^)') , (4.9) 
^ , 1 ^A-4/JN 27V2A^+41^^ 



J' ^ 4 \J'J 128 

A second order expansion of the dispersion relation at A = 1 has already been presented 
before [8], though with floating-point coefficients. Eq. (4.8) agrees with the result of [8] for 
the gap uo~{k = n) up to first order, but there is a minor difference in the second order: We 
believe that the coefficient of cos2/c in eq. (23) of [8] should read —0.52781 . . . (not —0.469). 
We have also checked eq. (24) loc.cit. and in this case found perfect agreement. 

Given the low order of the series (4.8-4.10) we do not try to draw conclusions for the 
weak-coupling region from them. We have restricted to only second order since a sym- 
bolic computation of higher orders is very difficult. This is due to the many square roots 
encountered, as is already indicated by the results presented here. 



V. THE STRONG-COUPLING EFFECTIVE HAMILTONIAN OF A 

FRUSTRATED LADDER 



Here we look at strong couphng (J' » J) for PBC and odd A^. In this case additional 
degeneracies prcchidc a simple analysis as in the preceding section. From a first-order 
consideration in J one infers that the low-energy effective Hamiltonian for (3.1) with A = 1 
and /i = is then given by (see [54,52] for A^ = 3 and [57] for larger A"): 
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HeS.'^^ = ly + «^ {(^t(^x+l + f^x+l)) S^S^+i , (5.1) 

x=l 

where the Sx are su{2) operators acting in the spin-space and act on another two- 
dimensional space which comes from a degeneracy due to the permutational symmetry of 
the chains. 

We have checked the vahdity of (5.1) for = 3, 5, 7 and 9 in the foUowing way: First 
one has to determine the groimdstatc space at each rung for J = which is nothing but the 
groundstate space of an A'^-site Heisenberg chain. For N odd, the lowest energy states have 

— ±1/2. In the case of OBC, this would be the only degeneracy. SU{2) symmetry is 
then sufficient to conclude that the effective Hamiltonian is a simple Heisenberg chain which 
is gaplcss in accordance with the generalized Haldane conjecture. 

For PBC there is another two-fold degeneracy in addition to this two-fold degeneracy 
in spin space: For N odd and PBC the groundstates of a Heisenberg chain carry momenta 
A; = ±27rL(iV + l)/4j/A^ where parity symmetry is reflected in the freedom of choice of sign. 

— * 

So, the groundstate-space at each rung x is 4-dimensional: The operators Sx act in the 
two-dimensional spin space and the cr^ act in the two-dimensional space spanned by the 
groundstate momenta. 

This degeneracy makes perturbation expansions in J highly non-trivial: At first order 
in J one has to diagonalize the matrix (5.1) which is determined by the matrix elements of 
the interaction along the legs in (3.1). That the only non-zero matrix elements are those 
given in (5.1) can be inferred just from the following symmetries of the full Hamiltonian: 
Global SU{2) symmetry (actually one needs only the U{1) Cartan subalgebra of su{2)) 
and invariance under simultaneous translations or refiections along all the rungs. These 
symmetries also imply some identities between the non-zero matrix elements, but at the end 
one still has to explicitly compute some matrix elements - at least in order to determine the 
constants ctjv- We have performed such direct computations of matrix elements for N — 3, 
5, 7 and 9 and found the associated values of ccjv to be 

^3 = 1 , "5 = ^ , ^7 = 2.6206859 . . . , ctg = 3.5012083 .... (5.2) 
9 

In contrast e.g. to the XXZ-chain (2.1), already for A^ = 3 the Hamiltonian (5.1) does not 
satisfy the Reshetikhin criterion (eq. (3.20) on p. 101 of [58]). Therefore, it is in general 

not integrable (in the sense that it would be the Hamiltonian of a one-parameter family of 
transfer matrices which commute among themselves and with this Hamiltonian). So, one 
has to treat it by other approximate or numerical methods; e.g. a DMRG study was carried 
out for A^ = 3 in [57] providing evidence for a gap to S — 1 excitations. 

In the present paper we are interested in generic XXZ-anisotropies A 7^ 1 in the inter- 
action along the chains in (3.1) and thus we should generalize (5.1). This generalization is 
obvious from the way how the XXZ-anisotropy appears in (3.1) and our derivation of (5.1): 
A just multiplies the matrix elements of the S^-S^ interaction. Therefore, the effective 
Hamiltonian for generic A is given by 

H^a'"^ = ^ E (1 + «iv + {as:S^,^, + i {S^S-^, + 5-5++ J) , (5.3) 
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where the parameters aN remain those given in (5.2). The generahzation (5.3) includes in 
particular the case A = 0, corresponding to two coupled XY-modcls. Then {i.e. for A = 0) 
one can apply a Jordan- Wigner transformation to (5.3). However, even in this case one 
obtains a four-fermion interaction with the effect that the problem does not simplify (in 
contrast to the familiar case of fermion bilinears). In particular, the determination of the 
groundstate of (5.3) for A = is far from being straightforward. 

As was pointed out in [15], the effective Hamiltonian (5.3) describes the response of (3.1) 
to a magnetic field for |(M)| < 1/N at strong coupling. For N = 3 {i.e. = 1) and A = 1 
we find by exact diagonalization of (5.1) that the transition to (M) = 1/3 (full magnetization 
for the effective Hamiltonian) takes place at 3h/J = 4.3146, 4.3121, 4.3108, 4.3100, 4.3096 for 
L = 8, 10, 12, 14, 16, respectively. This is in reasonable agreement with numerical values for 
the lower boundary of the (M) = 1/3 plateau of (3.1) at J' ^ J (compare Fig. 4 of [15]). 

VI. NUMERICAL ANALYSIS OF THE STRONG-COUPLING EFFECTIVE 

HAMILTONIAN 

To learn more about the spectrum of (5.3), we have performed numerical diagonalizations 
mainly for = 3 on finite systems, as was already done in [57] for A = 1. The Hamiltonian 
has two conserved quantities: (for A = 1 actually the total spin S) and a second similar 
quantity related to the first factor in (5.3) which wc denote by S^. The lowest eigenvalues 
are located in the = sector. First we look at the gap to the excitations in the = 1 
sector. It turns out that one can fit the system-size dependence of this gap nicely by ^ 

Es^=o,s^=i(-^) = Ee^=o,5^=i(oo) + ^ . (6.1) 

Estimates for these parameters based on data for lengths up to L = 14 are presented in Table 
HI for some values of A and N — 3. The numbers in brackets indicate the Icr-confidence 
interval of the fit for the last given digit. Since this ignores possible other error sources, 
the true error may be a little larger. Our result for A = 1 (£2^=0,5^=1(00) = 0.208(1) J) 
agrees within error bounds with that of [57] (£2^=0,5^=1(00) = 0.27(7) J). From Table HI 
we conclude that the = 1 excitation of (5.3) is gapped for all A > and that there is 
nothing special about the case A = 1 from this point of view. 

One comment is in place regarding the form (6.1) since in a gapped situation the conver- 
gence should ultimately be exponential (or at least of order - see e.g. [55] and references 
therein). Here we seem to observe a typical crossover phenomenon, i.e. the small values 
of the gap imply a large correlation length such that our system sizes may be well below 
the correlation length. In such a range of system sizes one would indeed expect to observe 
finite-size corrections which are typical for massless situations. Since the corrections should 



^This appears to be also a good fit to the data presented in Fig. 1 of [57] at least for intermediate 
system sizes and would lead to results compatible with the exponential fit used there. Note also 
that convergence is typically faster and clearer for the periodic boundary conditions used here than 
for open boundary conditions (used in [57]). 
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ultimately become smaller, this would lead to obtaining systematically too small values of 
the gap. With the fit (6.1) we thus obtain a lower bound for the gap which is presumably 
not far from the true value. In particular, we can safely infer the presence of a gap. 



A 


0.2 


0.4 


0.8 


1.0 1.2 


1.0 


Es^=o,5^=i(oo)/J 


0.139(5) 0.134(8) 


0.166(5) 


0.200(1) 


0.208(1) 0.214(2) 


0.390(6) 




2.78(1) 3.35(6) 


3.19(4) 


3.079(7) 


3.088(7) 3.13(2) 


5.72(1) 



TABLE III. Parameters for the fit (6.1) to the = 1 gap of (5.3). The first six columns are 
for N = 3, but various values of A. The rightmost column is for 'W = cxd" with un/N —>■ 1. 



Concerning the case of N > 3, one observes from (5.2) and a further value for an [57] 
that aN is roughly proportional to N, i.e. ^ 0.44 A?" for large N. Using this information, 
we have extrapolated (5.1) to infinite N setting MmN^ooCtN/N = 1 in order to avoid the 
uncertainty in the true proportionality constant. This limit eliminates the term 1 x SxS^+i 
in (5.1). The rightmost column in Table III shows the value for the S* = 1 gap that we obtain 
in this case. It should be noted that the properly rescaled value for iV — cxd is slightly lower 
than that for = 3 at A = 1 (the former is about 80% of the latter). However, even for 
N — oo our estimate for the gap is still remarkably distinct from zero. This suggests a gap 
in the strong-coupling hmit (5.1) for all N which slightly decreases as N ^ oo, but does not 
close even in this limit. 

Now we turn to the "gap" in the = sector for N = 3. The data in Table IV can be 
interpreted as evidence that it asymptotically decreases roughly as 

Ee^=o,5^=o ~ j2 ; (6-2) 

at least this "gap" clearly tends to zero in the thermodynamic limit (in particular close to 
A = the finite-size exponent could be different from that given in (6.2)). This energy 
level corresponds to the state constructed in the generalized Licb-Schultz-Mattis theorem 
[21-23,12]. According to [57] this energy- level should be interpreted as a degenerate ground- 
state arising due to spontaneous dimerization in the thermodynamic limit. The fact that 
the energy levels in Table IV have momentum tt relative to the groundstate is compatible 
with this interpretation and yields Z = 2 for the condition (1.1). 
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11.04294 


12.16715 


13.29556 


15.54226 


16.65656 


17.76262 
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14 
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21.60791 


23.19367 



TABLE IV. Rescaled 5^ = "gaps" of (5.3) with iV = 3 for various values of A. 



Next, we investigate the momentum-dependence of the gaps to the lowest excited states 
of (5.1) with N — 3. The data for = and total spin 5" = is shown in Fig. 3 and that 
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for total spin S = 1 (also = 0) in Fig. 4 (compare also Fig. 4 of [57]). Here, wc measure 
the momentum of the excitations relative to the groundstate. It should be noted that due 
to parity conservation only half of the spectrum is shown (the parts for A; > tt or A; < are 
mirror-symmetric extensions of this figure). The two figures look quite similar. Both can be 
interpreted as the lower boundary of a (two-particle) scattering continuum, In particular, 
we do not observe one-particle states. 
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FIG. 3. Lowest gaps of (5.1) with = 3 in the sector with = and total spin S" = as a 
function of momentum k relative to the groundstate. The symbols are for L = 6 (O), L = 8 (+), 
L = 10 (□), L = 12 (x) and L = 14 (A), respectively. The line is the extrapolation (6.3) of the 
lower boundary L ^ oo. 

To extrapolate the lower boundaries of these two-particle scattering states, we have 
Fourier transformed ^. Then we have extrapolated each coefficient of the Fourier series 
separately using a Shanks transform (which is the a — special case of the vanden Broeck- 
Schwartz algorithm - see e.g. [59]). This leads to: 

^l.=o,s=oik) / J = 0.654(478) - 0.014(191) cos A; - 0.411(108) cos2A; 

+0.040(136) C0S3A; - 0.044 cos 4A; , (6.3) 



^We Fourier transform rather than E, since then higher harmonics are better suppressed, as 
must be the case for the scaling limit of a lattice model close to a second order critical point. 
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=o,s=iik)/J = 0.671(43) + 0.023(24) cos A; - 0.492 cos 2A; 

-0.025 cos 3k - 0.058(11) cos 4A; - 0.013(50) cos 5A; 



(6.4) 



The numbers in brackets indicate estimates for the error of the last given digits. Here we 
have suppressed the highest harmonics, since they cannot be reasonably extrapolated, but 
are expected to be small anyway. 

In this way we obtain a rather inaccurate estimate for the gap E ^ 0.3J with a large 
uncertainty which is due to the large errors in particular in (6.3) and the uncertainty in 
the higher harmonics. Nevertheless, this estimate is still quite close to the one in Table III. 
A more interesting observation is that (6.3) and (6.4) are equal within error bounds. This 
suggests that these two thresholds can be interpreted in terms of two-particle scattering 
states of a single fundamental particle. Such a fundamental excitation would have to be 
similar to the spinon in the XXZ-chain; in particular it would have to carry S = 1/2 (and 

= ±1/2). 




FIG. 4. Same as Fig. 3, but for total spin S = 1. The line shows the extrapolation (6.4). 



Let us now try to exhibit this fundamental excitation explicitly. For even L and peri- 
odic boundary conditions we have only found two-particle scattering states in the low-lying 
excitation spectrum. Therefore, it is natural to look for a spinon-type excitation at odd L 
(still periodic boundary conditions) in the same way as one can exhibit the spinon for the 
XXZ-chain [38-40]. We have computed the spectrum of (5.1) for N = 3 and odd L from 
5 to 13 in the sector with = 1/2 and total spin S — 1/2. The main difference between 
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the present situation and the XXZ-chain is that here we expect a charge conjugate pair of 
spinous (S^ = ±1/2) while for the XXZ-chain there was only one. 

Making single-particle states visible is traded for the absence of a groundstate at odd L. 
In order to permit interpretation of the results as gaps we have therefore interpolated the 
groundstate energy using the values at L± 1. The resulting dispersion relation for the spinon 
is shown in Fig. 5. As was already the case for the spectra at even L, it turns out that k 
should be defined such that translationally invariant states on the lattice appear alternately 
at A; = and k = tt - the actual convention can be read off from Fig. 5 noting that only 
either A; = or A; = tt can be realized for odd L. 

1.8 n 1 1 1 n 

o 

1.6 - 




u i i i u 

71/4 71/2 3 71/4 71 

k 

FIG. 5. The spinon of (5.1) with = 3, i.e. lowest gaps in the sector with = 1/2 and total 
spin S = 1/2. The symbols are for L = 5 (<>), L = 7 (+), L = 9 (□), L = 11 (x) and L = 13 
(A), respectively. The line is the extrapolation (6.6) of the dispersion curve to the thermodynamic 
limit. 

To interpret the data, we have again Fourier transformed E^. Firstly, this gives an 
interpolation of £2^=1/2,5=1/2 at /c = 7i/2. Analogously to (6.1) we fit the data for L — 5,9, 13 
to the form (the values for L = 7 and 11 should be omitted to obtain a monotonic sequence) 

Es^=l/2,5=l/2(-f') = Esz=i/2,5=i/2(oo) + — , (6.5) 

and obtain an estimate for the gap of the spinon £2^=1/2,5=1/2(00) = 0.131(8) J with a — 
0.51(6) J. This is roughly consistent with half the value in Table III or that given in [57], as 
it should be if our interpretation as single- respective two-spinon scattering states is correct. 
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An alternate way to analyze the data is to extrapolate each coefficient of the Fourier 
series separately using a Shanks transform. Using now all available L, we find 

^h=i/2,s=i/2ik) / J = 0.6331(155) + 0.0592(184) cos A; 

+0.5387(122) cos 2k + 0.0121(63) cos 3k (6.6) 
-0.0633(160) cos4:k - 0.0127 cos 5A; + 0.0177 cos 6A; . 

As before, the numbers in brackets indicate estimates for the error of the last given digits. 
For the two highest harmonics there is not sufficient data for an extrapolation, so we just 
take the L = 13 estimate without being able to estimate an error. The extrapolation 
(6.6) is shown by the line in Fig. 5. Obviously, finite-size effects are more important for 
k < 7r/2 than for k > 7r/2. Eq. (6.6) yields another estimate for the gap of the spinon 
Es^=i/2,s=i/2(oo) ~ 0.116 J. The error estimate obtained from (6.6) is not sensible, but the 
value for the gap itself is very close to our previous extrapolation or half the value given in 
Table III. 

Finally, we have checked that within error bounds the dispersion relation (6.4) can be 
written in terms of (6.6) as Esz=o,s=i(^) = £22=1/2,5=1/2(^ — ^0 + Es^=i/2,s=i/2(^"') with some 
k'. Such a decomposition must be possible if our particle interpretation is correct. 

Methods similar to the ones used in the present section may be useful also in other 
cases beyond the present one and the study of [40]. One natural such candidate is a direct 
observation of a spinon-type excitation in = 3 cylindrically coupled chains at intermediate 
or small couplings J'. 

VII. SUMMARY OF RESULTS 

Our results are best summarized in (schematic) magnetic phase diagrams. For definite- 
ness we consider the SU{2) symmetric situation A = 1, though similar pictures can be 
drawn for other values of A as well. 

For completeness, let us start with the case N = 2, where the corresponding picture is 
given by Fig. 6a). The boundary of the (M) = plateau is determined by the spin-gap in 
zero field which for an A^ = 2 leg ladder has been studied in great detail. For J'/ J < 1.5 we 
use the Quantum Monte Carlo results of [11] in Fig. 6a); for J' / J > 1.5 the raw 13th order 
strong-coupling series of [56] is used instead (note the excellent matching at J' / J — 1.5). 

Both the numerical data [11] and the series expansions [56] support a linear opening of 
the gap for small J', as was predicted by a dimensional analysis of the perturbing operator 
in the field-theoretic formulation [60-62]. 

Fig. 6b) shows a more interesting case, i.e. N — 3 with OBC. The boundaries of the 
(M) = 1/3 plateau have been determined from the fourth order strong-coupling series (4.3) 
and (4.4) for J'/ J > 2. For 1 < J' / J < 2 we obtained them from a Shanks extrapolation 
of the finite-size data in our earlier paper [15]. The remaining weak-coupling region is 
the most speculative part of the figure. We located the ending-point of the (M) = 1/3 
plateau in the vicinity of the corresponding point on the magnetization curves of decoupled 
Heisenberg chains, as is suggested by the Abelian bosonization analysis if one assumes a 
similar behaviour for OBC and PBC (the bosonization analysis predicts for PBC that the 
(M) = 1/3 plateau disappears for small but non-zero J'/J)- 
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FIG. 6. Schematic magnetic phase diagram at A = 1 for a) iV = 2, b) A?^ = 3 and OBC, c) 
= 3 and PBC, d) = 4 and OBC, e) AT = 4 and PBC. White regions in the h-J' j J plane 

indicate gapped regions with a plateau in the magnetization curve while the shaded areas are 
massless and the magnetization (M) changes continuously if the applied field h is varied. 
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The analogous case with changed boundary conditions {i.e. N — 3 and PBC) is shown in 
Fig. 6c). Here, no series expansions are possible due to extra degeneracies at strong couphng. 
The boundaries of the plateaux have therefore been determined in this case mainly on the 
basis of older numerical data [15]. We have used a Shanks extrapolation for L = 4, 6 and 
8 for J'/ J > 2.75 at the lower boundary of the (M) — 1/3 plateau and for J' / J — 2.5 at 
its upper boundary to estimate their location. For smaller couplings, the finite-size data is 
non-monotonic. The best we can do in the range 1 < J'/ J < 2.5 is to fit the L = A and 8 
data to a form with 1/L corrections like (6.1,6.5). The ending point of the plateau is again 
placed on the basis of the weak-coupling analysis. It should be noted that the nature of 
the transition at the upper boundary of the (M) = 1/3 plateau changes qualitatively for 
J'/ J > 2.75 [15]. This strongly frustrated region is indicated by the bold line in Fig. 6c). 
It is possible that the transition becomes first order along this line. Note that the region 
in question is far outside the weak-coupling region, where we expect all transitions to be 
continuous. 

Another interesting difference between Fig. 6b) and Fig. 6c) is that in the latter a tiny 
(M) = plateau {i.e. a gap) opens. Its boundary has been estimated at intermediate 
couphngs by fitting the L = 4, 6 and 8 data [15] to the form (6.1). This yields shghtly 
smaller values than those given in [57] in the cases where we overlap. However, we agree 
with [57] in the most important point, namely the existence of such a plateau. A numerical 
determination of its ending point is difficult, and the field-theoretical weak-coupling analysis 
is not yet conclusive either [53,54] - the ending point may well be anywhere in the weak- 
coupling region < J' < J. 

Finally, the magnetic phase diagrams for = 4 are given by Figs. 6d) and e) for OBC 
and PBC, respectively. To obtain them, we have performed further numerical computations. 
For the upper boundary of the (M) = 1/2 plateau we have numerical data for L = 4, 6 and 
8 such that we can apply a Shanks transform to it. For its lower boundary we have only 
L — 4 and 6 data and therefore we have to make an assumption on the finite-size corrections 
to extrapolate it (we assumed 1/L corrections, though this is not entirely satisfactory). In 
the strong-coupling region we used our series instead of numerical data. The series and 
numerical data are matched at J' / J — 3.5 or J'/ J — 2.5 for the series corresponding to 
the upper boundary for OBC (4.10) or PBC (4.7), respectively. At the lower boundary 
of the (M) = 1/2 plateau we matched the series (4.9) and (4.6) to the numerical data at 
J' /J = 2.25 and J' /J = 1.75, respectively. Neither of the methods accessible to us is very 
accurate in the region where this plateau closes, but all three methods (numerical, series 
and Abelian bosonization) point to a location of the ending point in the region where J' and 
J are of the same order. 

The gaps for A^ = 4 are taken from our series (4.8) for J'/ J > 2 for OBC and (4.5) 
for J' / J > 1.5 for PBC. For OBC the accurate numerical data for the gap of [11] is used 
in the weak-coupling region. The corresponding line in Fig. 6e) is in comparison rather an 
educated guess which is inspired though by L = 4 and L — 6 numerical data. Regarding the 
series both for the gap and the boundaries of the (M) — 1/2 plateau, we observe a trend 
that those for PBC can be used for somewhat smaller values of J'/ J than those for OBC. 
This is expected since the former are fourth order but the latter only second order. 

Although Fig. 6 is for the particular choice A = 1 there is nothing particular about this 
case (at least for non-zero magnetizations), and one would obtain similar figures for other 
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values of A as well. Further plateaux may open for A > 1. In particular, there should 
always be an (M) = plateau in N coupled XXZ-chains with A > 1, since each such chain 
is massive and this should be preserved at least for sufficiently weak coupling. In the Ising 
limit A — > oo and for non-frustrating boundary conditions it is easy to see that this is 
accompanied by breaking of translational symmetry to a period I — 2 in the groundstate. 
In the general case A > 1, such a period I — 2 reconciles the appearance of a gap for both 
even and odd N with (1.1). 

The Abelian bosonization analysis predicts all the massless shaded regions in Fig. 6 to 
be c = 1 theories (with the exception J' — where one trivially has a, c — N theory). In 
these regions the exponents governing the asymptotics of the correlation functions depend 
continuously on the parameters. Predictions can be made, however, for the transitions at 
the boundaries between such massless phases and plateau regions. The opening of a plateau 
when varying J' is a transition of K-T type [45] . Like in the case of the transition at A = 1 
in the XXZ-chain, this implies a very narrow plateau after the transition (c.f. (2.12)) which 
makes it difficult to observe numerically [15]. At the transition point the asymptotics of 
the correlation functions is governed by the exponents (3.14) while along the boundaries of 
the plateaux one has the universal exponents (3.15). It should be noted that an attempt 
to verify the latter exponents numerically or experimentally is likely to rather lead to the 
exponents characteristic for the transition point if one is sufficiently close to it. 

The field-theoretical analysis also predicts the asymptotic behaviour of the magnetization 
in a massless phase but close to a plateau boundary to be given by the universal DN-PT 
behaviour [42,43] (2.10). We have in fact numerically verified such a square-root behaviour 
close to saturation ((M) — ^ 1) at some values of J' /J for iV = 2, 3 and 4 with both 
OBC and PBC. However, the example of the XXZ-chain shows (c.f. Fig. 2) that close to 
other plateau boundaries this universal behaviour may be restricted to a tiny region and its 
observation could be very difficult. In experimental situations, it will be further obscured 
by thermal fluctuations and other effects such as disorder (see e.g. [17]). This explains why 
rather accurate experiments on = 2 leg spin-ladder materials [16,63] show no evidence of 
a square-root behaviour for (M) — > 0. 

Quite surprisingly, massless excitations (though non-magnetic ones) also arise in plateau 
regions. This can be seen from (5.3) which for (M) = 1/A^ is just an XY-chain and therefore 
massless. This yields massless excitations in the limit J' ^ oo in Fig. 6c), or more generally 
in the strong-coupling limit on the (M) = 1/N plateau for A^ odd and PBC. Whether 
J/J' = is just a critical point or if massless non-magnetic excitations also arise at finite J' 
remains to be investigated. 

VIII. DISCUSSION AND CONCLUSION 

In this paper we have investigated the conditions under which plateaux appear in A^-leg 
spin ladders as well as the universality classes of the transitions at the boundaries of such 
plateaux. Certain small plateaux may have slipped our attention. For example, there could 
be a narrow (M) =2/3 plateau for A" = 3 and PBC at intermediate or strong couphng 
J' which would be accompanied by spontaneous breaking of translational symmetry to a 
period 1 — 2. If this should turn out to be the case, it would have to be added to Fig. 6c). 
However, our main point is the presence of such plateaux, not the absence of particular ones. 
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We also confirmed the conclusion of [57] that in the case N — 3 frustration induces a 
zero-field gap at least for sufficiently strong coupling. It may be even more intriguing that, 
according to our strong-coupling data, this gap seems to survive the N —>■ oo limit for an odd 
number A'^ of cylindrically coupled chains. This shows that it is necessary to specify at least 
boundary conditions along the rungs in the generalization of the Haldane conjecture to spin 
ladders. The peculiar behaviour of a cylindrical configuration may be interpreted as follows: 
Strongly frustrating boundary conditions force a one-dimensional domain wall into the two- 
dimensional system corresponding to iV = oo. The N ^ oo limit of the one-dimensional 
Hamiltonian (5.3) is just the effective Hamiltonian for the low-energy excitations of this 
domain wall. As a consequence, there cannot be any long-range order which is typical for 
the two-dimensional Heisenberg model, and there is no reason why the low-energy spectrum 
should not be gapped as it apparently is. 

Similar surprises cannot be completely ruled out in the weak-coupling region for even 
and (M) = 0. The zero-field case is difficult to control since there is an additional relevant 
interaction between the massive degrees of freedom and the possibly massless ones (the 
coefficient of A2 in (3.4)). In general, this gives rise to non-perturbative renormalization. In 
the isotropic case A = 1 one can use the SU (2) symmetry [52] to infer the renormalized 
radius of compactification of the remaining massless field. This then leads to the generalized 
Haldane conjecture in the framework of Abelian bosonization. For A > 1 a gap is always 
expected in the weak-coupling regime since already the decoupled chains are massive. The 
situation is far less clear for anisotropics A < 1 since then it is not known how to control 
the renormalization group fiow. An investigation of this region by other non-perturbative 
methods would be interesting. It may also be desirable to perform further checks of the 
absence of an extended massless phase in the SU{2) symmetric situation for even A^ > 6. 

Beyond a more detailed understanding of the A^-leg spin-ladder model (3.1) treated here, 
a similar investigation of other models could be interesting. One natural step would be to 
include charge degrees of freedom, and see if interesting effects arise from the interplay of a 
magnetic field with transport properties. 

Last but not least, it would be desirable to have an experimental verification of our 
predictions. We are confident that this is possible in principle and hope that it will in fact 
be carried out. 
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